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ABSTRACT 


Three different approaches to the problem of imple- 
menting a reduced-order, sub-optimal Kalman filter for a 
discrete, linear stochastic process, with time-invariant 
dynamics, are presented. 

A first method, A, is based upon the partitioning of 
the system dynamics. A second method, B, is implemented 
using matrix pseudo-inversion and a third method, C, is 
based upon reduction of the original process to one of 
lower order using the dominant roots of the system. An 
expression for the performance degradation in method A is 
derived. In method B, expressions for the sub-optimal 
estimation error, and sub-optimal variance of estimation 
error are derived. 

The several methods are applied to a fourth-order 


process for illustration. 
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I. INTRODUCTION 


In any modern control problem involving filtering of 
discrete or continuous stochastic linear’ processes, it is 
clear that the computational requirements associated with 
the implementation of a Kalman filter increase consider- 
ably as the dimension of the system increases. A simpler 
(reduced) configuration of the filter with sub-optimal 
performance may be acceptable, mainly if some preselected 
r states (where r=n, and n represents the order of the 
original system) are of primary interest in the estimation 
process and if the performance degradation incurred when 
using the reduced model is not severe. 

With this idea of implementing a reduced, discrete, 
sub-optimal filter, based upon the optimal Kalman estima- 
t#on theory (Ré&.).), a Dibfary séarch was done, yielding 
two papers of particular interest which present different 
approaches to the problem of reduced estimation for con- 
tinuous dynamic processes. Meditch (Ref.3) in 1964, de- 
veloped a theory of sub-optimal filtering for continuous 
dynamic processes, in which the process vector is divided 
by means of matrix pseudo inversion, with separate filters 
designed for each part. The sub-optimal estimates are ob- 
tained by combining the several filter outputs. As 
Meditch points out in his paper, an important feature of 
his continuous formulation is that it provides a powerful 


tool for preliminary analysis of sub-optimal filters. 


Huddle and Wismer (Ref.4) in 1968 based a sub-optimal 
filter design on a "primary system", which is implemented 
by dividing the original process into two separate systems, 
using partitioning of the continuous dynamics of the sys- 
tem, and neglecting the remaining “secondary system’. 

Another successful technigue for simplifying filter 
computations has been investigated by Aoki and Huddle 
(Ref.5), wherein a procedure is developed to estimate the 
state vector of a discrete stochastic system, when con- 
straints are imposed on the number of memory elements of 
the estimator. 

This thesis presents both an extension of two of the 
methods referred to above (Ref.3 and Ref.4) to the diserete 
case, and a third approach to the same problem based upon 
the use of the dominant roots of the full process. 

The study of the three different approaches to the 
reduced sub-optimal filter problem is presented in a co- 
herent form as follows. Chapter II describes the general 
discrete linear stochastic process. In Chapter III, the 
work of Huddle and Wismer (Ref.4) is extended to the dis- 
crete case, and the performance degradation of the re- 
duced filter is obtained. An example using a fourth- 
order process is solved estimating only two states of the 
full system, in order to illustrate the method. In Chap- 
ter IV, the Meditch (Ref.3) reduced filter is adapted to 
the discrete case, and recursive equations for the sub- 


optimal estimation error and sub-optimal covariance of 


estimation error are derived and used as a means of com- 
paring the sub-optimal with optimal performances. The 
same fourth-order system used in Chapter III, is used here 
as an example, and again, only two states are estimated. 
In Chapter V, another reduced filter is implemented, using 
the dominant roots of the complete system in order to pro- 
duce the dynamics of the reduced filter. Once more, the 
technigue is illustrated using the original plant of the 
previous examples. 

In the three examples presented as illustrations, the 
same complete process is used in order to compare the re- 
sults of the three different methods. In each one of the 
examples, a Monte Carlo simulation of the full process was 
performed, and optimal Kalman estimation of the whole 
state vector was accomplished in order to establish the 
performance degradation of the reduced filters. 

Appendix A presents the basic principles of matrix 
pseudo-inversion; more complete treatments of this sub- 
ject appear in Penrose (Ref.8) and Werther (Ref.9). 

Appendix B and C include the computer programs used 
in solving the examples. These programs were not intended 
to be computationally efficient either in time, or in 


storage. 








IL. THE DISCRETE (DISGRE REZ ED main E AR 
SROCHASTIG PRO@ESS WIT 
TIME-INVARIANT DYNAMICS 


A linear, continuous stochastic dynamic process may 
be described in general (Ref.7) by the vector matrix dif- 


ferential equation 


k(t) = A(t)x(t) + B(t) u(t) + w(t) (2.1) 
with output 
y(t) = H(t)x(t) | Cm 


where: x(t) is a (n x 1) vector, representing the state 
variables of the system 

u(t) is a (m x 1) vector (mn), representing the 
deterministic input to the system 

w(t) is a (m x 1) vector representing the stochas- 
tic input to the system 

A(t)eis a (n x8n) systemmmatrix 

BCE) is a GaP im)” dist rebut fon marrrx 

y(t) is a (p x 1) vector of system outputs 


H(t) is a (p x n) measurement matrix. 


It will be assumed that measurement of the output 
can be achieved only with the addition of noise, so the 


measurement equation is defined as 
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ie) = y(t) + v(t), 2.3) 


where v(t) is a (p x 1) vector representing the stochastic 
measurement noise and z(t) isc De X 1) vector of measured 
outputs. 

The process is visualized better by means of Fig.2.1, 
where the double lines indicate a signal flow character- 
istic of all matrix block diagrams. 

It'- is assumed for simplicity that the process is sta- 
tionary (time-invariant); that is, all the elements of 
BaewA(t), Blt) and H(t) matrices are constants. It is 
further assumed that the process is to be sampled at dis- 
crete times and that u(t) and w(t) are constant vectors 
during each sampling period. The general solution for 
Eq. (2.1) under these conditions is expressed as 


4 ' > t : 
x(t) = CME x(t.) sig OE € B[u +t w aT (2 ois) 


If the sampling time is assumed to be T seconds, Eq. 


(2.4) may be rewritten as 
T a ; 
x(T) = €°'x(t,) + ef cATaT Bu 4+ w] (2.5) 
O ne 
Defining 


Q(t) Sct , (2. 6a) 


IAL 








Fig e2. 1 ocueBihock daagram for a Gonkandousmstoehast ie 


process with time-varying dynamics. 





Fig. 2.2 - Block diagram for a discrete stochastic 


process with time-invariant dynamics. 
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and 


P (Dee (2.6b) 
where O0(T) (n x n) and I'(T) (n x m) are called the dis- 

Sreveuctate Eransition and distribution matrices, respect- 
ively. The Picea miDeae EOEnM Of sd. (2.1) becomes in the 


stationary case 
x(k+1) = O(T)x(k) + I(T) (ul) + w(k) | ; Cn) 


0(T) and I(T) may be obtained by application of the 


Laplace transformation, i.e., 


O(T) -£7 (31 ee (2.8a) 


Veer 


Distinction is made between a truly discrete system,and 
one which is continuously operating but driven and measured 
at discrete times. The linear difference equation for the 
truly discrete system will be 


x(kt1) = Anx(k) + Bp u(k) + wk) 


The terminology "discrete", ‘discretized" is used without 
Gistinction from this point on. 
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ia.) = p+ O(s)B : . (2.8b) 


The discretized, linear stochastic process with time- 
invariant coefficients will hereafter be represented by 


the model (see Fig. 2.2) 
x(k+1) = Ox(k) + Dfu(k) + wk] (2.9) 
with measurement equation 
2 Cke=eHx (Ck) tek)» . (2.16) 


Eqs. (2.9) and (2.10) will be used throughout this 
thesis as the model of the complete discrete system for 
which reduced sub-optimal filters are implemented. 
Furthermore, w(k) and v(k) will be considered to be vec- 
tors of independent, gaussian random processes (white 
noise) with zero means and fully described by the follow- 
ime statistits: 


Elw(k)| = 0 for all k , (2.11) 


) 
t — 


1 


E.w(k)+w'(i)] = Efw(k)| -E[w'(i)]) = 0 — k # i, (2.12) 


Hi 


known’ quantity for all k , (2.13) 


Elw(k)+w'(k)) 
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E |v (i)! a) 
E[w(k)-v°(i)} = 


and 


E[v(k)-w'(k)} = 


where R is a (p x p) 


noise. 


formal ik , £4) 


E[v(k)] -E[y'(i)] = 0 k # i, (2.15) 


R fom all k . OP Us) 


matrix of covariance of measurement 


1) 


IIL. IMPLEMENTATION OF SA @RERUCED DISGEETES 
LINEAR KALMAN FILTER BY PARTITIONING 


THE DYNAMICS OF THE @RIGINAL PROCESS 


This chapter presents a discrete solution to the prob- 
lem of implementing a sub-optimal filter by means of par- 
titioning the original discrete system (Fig. 2.2) into two 
reduced models, a "primary system" model which will include 
only a selected subset of the state vector components whose 
estimation is of primary interest, and a "secondary sys- 
tem which will be formed by the remaining state varibales. 
Then, the primary plant will have reduced dimensions and 
the implementation of a reduced Kalman filter based upon 


this primary model will be possible. 


A. DERIVATION OF THE REDUCED FILTER EQUATIONS 

It 1s assumed that the system dynamics are represent- 
ed by Eg. (2.9) and that the statistics of the forcing vam 
measurement noise are fully described by Eqs. (2.11) 
through (2.16). For simplicity it will be further assumed 
that there is no deterministic forcing function u(k). The 


stochastic process may then be expressed as 
x(kt+1) = Ox(k) + Dwlk) , (300) 
with measurement 


z(k) = Ext awn 


lo 


By matrix partitioning, the discrete dynamics of the 


full process may be rewritten as 





(322) 


hephenemEmll state vecebor peine an (n x 1) matrix, 
then: 

X, represents an (ny x 1.) state vector which includes 
all those states of particular interest in the estimation 
process; 

Xy represents an (nj x 1) vector of the remaining 
state vector components (ny ele n); 

b, 1s an (n, x n,) matrix, representing the dynamics 
of x3 

0. is an (ny x ny) matrix, representing the dynamics 
Or X93 

0.4 and 0.9 are (ny x no) and (n, x ny) matrices re- 
spectively, indicating the coupling between the primary 
and secondary systems; 

w, is a (m, x 1) vector of forcing (white noise) 
mime t ions; 

Wy is a (m, x 1) vector of forcing (white noise) 
functions with my) =r mM, = 0; and 

I) and I, are (ny Xx m,), (n. x mM, ) distribution ma- 


trices respectively. 


IE? 


The observation matrix H is partitioned as follows, 


The measurement equation may therefore be rewritten as 


2(k) = | Hy! Hy. P+ v(k) (3.3) 


It is assumed that all the measured outputs are in- 
cluded in the primary system, which makes H, = 0. This 
partitioning may not always be possible, and even if it 
is, it may require special rearrangement of the z vector. 

In Eq. (3.3) z(k) is the original (p x 1) measured 
vector and Hy isethe wGpe x ny) measurement matrix for the 
primary system. It may be seen from this special parti- 
tioning of the H matrix that the primary system vector 
xX] will not only be formed by those states of primary in- 
terest in the estimation process, but will also include 
those states whose estimation may not be needed, but 
whose outputs are being measured. This will establish 
the’ ComGd2tt One Waits n,=P- 

As a result of the preceeding discussion, Eq. (3.3) 


becomes 


z(k) = H)x,(k) + v(k) , (3.49 
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which represents the measurement equation for the reduced 


system. 


The general system may be described as consisting of 


a primary system. 


x (k+1) = O)x,(k) + O.yx5(k) + Tywy(k) , (3. 5a) 


and a secondary system described by, 


Xo(k+1) = O.5x,(k) + O)x5(k) + Thwo(k) . (3.51) 


Fig. (3.1) shows a complete representation of the parti- 
tioned system. 

In order to reduce computational requirements in the 
estimation of the preselected states, the design of the 
filter is based on the dynamics of the primary system, 
which can be represented as a simplified model of the re- 


duced process by 


x, (kt1) a 0, x(k) Be Tw, (k) ) : (3.6) 


1 
A bar over a quantity indicates a reduced-model quan- 
Peet y . 
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with measurement 
2(k) = Hyx,(k) + v(k) 


In the simplified model of the primary system the term in- 
volving b. 1X9 (k) has been neglected due to the fact that 
the secondary state vector X(k) is not considered at all 
in the implementation of the reduced filter. The conse- 
euences Of using Eq. @3.6) instead of Eq. @@@5a) will be 
analyzed in Section B of this chapter. 

If the Kalman filter algorithm is applied using the 
dynamics of the full plant, the matrix equation for the 
gain G(k), the covariance of estimation error P(k/k) and 
the prediction covariance P(k+1/k) are given in the dis- 


crete case by 


Gk) = P(k/k-1)H! 4H P(k/k-1)H* + R ; (3.7a) 
P(k/k) = I - G(k)H P(k/k-1) , (3.7b) 
and 
P(k+1/k) = @ P(k/k)O' + Q , (3.7c) 
where 
Q= VE(w(k)w(k)]P*. 


G(k) is an (m x p) matrix, P(k/k), P(k/k-1) and P(k+1/k) 


are (n x n)matrices. 


Jd 





The estimation equation is given by 


R(k/k) = O8(k-1/k-1) + 6(k)'2(k) = HOR(k=1/k-1). 


(3.39 


where X(k/k) is an (n x 1) vector of estimated states. 
By induction, the reduced sub-optimal expressions 


using the dynamics of Eq. (3.6) are 
Gal@k ) e=aar (k/k-1)H, [H P Ck lee WH + . « (3.9a) 
1 1 Lalani iL in eS 
: ; 
P, (k/k) nial oc G1 (k)H,. Py (k/k-1) ; (3.9b) 
and similarly 
IL 
P,(k+1/k) = 6,P,(k/k)d; + Q, , (3.9c) 
where 


Q, = TElwy(k)wi(k) rT 


P1(k/k) (ny x n,) represents the sub-optimal covariance 
of estimation error and P, (k+1/k) (ny Xx ny) is the sub- 


optimal prediction covariance of estimation error. 


ZZ 


The sub-optimal estimation equation is 


1 (k/k) = , x) (k-1/k-1) + G,(k)[2(k) - yO) x) (k-1/k-1)). 
(3.10) 


Considering the problem of initialization for the 
discrete reduced covariance equations, it is seen as a 
good approach to use for P,(0/-1), the partitioned por- 


tion of P(O/-1) corresponding to P,,(0/-1), or considering 


that 
Sy, : 
P,,(0/-1) | P,,(0/-1) 
P(0/-1) = |---~-------- 
egg e! ~1) | P9(0/-1) 
then 
he) = P77 (0/-1) . 
Hersthe filter initial states, 
= “A 
xj(0/-1) = x, (0/-1) , 
ix 
where 


*%(0/-1) 


&,(0/-1) 


ZS 





Fig. 3.2 shows a complete block diagram of the reduced 


plant and filter. 


S! REDUCED FILTER DEGRADATION 
Because the dynamics of a simplified system repre- 


sented by 
X,(k+1) = O)x, (k) + Tw, (k) (3.6) 


were used for the implementation of the reduced filter, 


instead of the dynamics of the primary system 


Sl) 0, x, (k) + 0. 1X9 Ck) + Tyw,(k) , (3.53) 


the estimation of xX, Is caused to be sub-optimal and in 
error. This error may be considered the performance de- 
gradation of the sub-optimal filter. 

Considering the reduced filter, its estimation error 


covariance is given by 


Vk) = E) xy (ke) = (k/e)] (3.11) 


where 


x, (k) = 0) x) (k-1) + 0.) x5(k-1) + Tywy(k-1) , (3. 5a) 


and by proper use of Eq. (3.6) it may be rewritten as 


x(k) = X)(k) + Oy X9(k-1) . (3.1 
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In the same way, 


x, (k/k) = §,) (k-W/k-1)+6,(k) (20k) - HO, x) (k-1/k-1) | 
3, 0) 
which by substitution of Eq. (3.4) for z(k), gives 


)(k/k) = Ox, (K-1/k-1)+ 6, (k) {Hy xk) +w(k) - Hy G,3(k-1/e) } 


iiegrhetssubstitutien sof He. (3282) for x, (k) gives 
x, (k/k) = §,x,(k-1/k-1) + 6, (Hy | () + 9. 4X (k-1) 


- 618 (k-1/k-1) | + G, (k)v(k) 
(C3eallag) 


Using Eqs. (3.12) and (3.13), the experimental estimation 
error covariance V(k) (Eq. 3.11) may be rewritten, after 
proper rearrangement, as 


V(k) = E 


(1, - 6G, (k)H, (x, Ck) - Q) x) (k-1/k-1) 


L 


+ (1, = GC) Hy 5, p39 (KL) + 6, (kdw(K)] "|, 
(3.14) 


By application of matrix algebra and elemental algebra of 
expectations, Eq. (3.14) may be expanded and after some 


manipulation gives 


Z> 
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V(k) = Ty -G) (K)H, | Py (k/k=1)+ 1, -G(k)Hy | 9.8 [x5 (k-1)| O01 


+8 [By (k/k-1) x9 (k-1)] O07 aE [xq DSfeeAe-1)] | fy -G, CDH, |" 7 
(3.15) 
where 
Ey (k/k-1L) = Ky (k) = Hy (k-1/k-1) 


The first seem on time wright hand side of Eq. (3.15) 
can be easily recognized as P,(k/k), the sub-optimal 
theoretical variance of estimation error based on the re- 
duced model. It is intuitively expected that P,(k/k) > 
P11 (k/k). Tiberaehiasim ne seins VOL wanes lia or epresentyethe 
further degradation of the reduced filter and are seen to 
be functions of that partitioned portion of ) which couples 
Xo to the primary system, i.e., bo: This suggests that as 
the coupling between x» and the primary system (O04) de- 
creases, the quantity P, (k/k) more accurately represents 
the actual estimation error covariance. 
oe EXAMPLE 1 

The process used for this example was a fourth-order 
model, with continuous linear stationary dynamics subject- 
ed to a forcing function having both a deterministic and a 
random component, with the latter having zero mean and 
known variance. It was assumed that only one state vari- 
able was measured, subject to the presence of gaussian 


measurement noise with zero mean and known variance. at 


27 


was further assumed that both random sequences have sta- 
tistics as described by Eas. (2.11) through (25!) —mae 
reduced filter was implemented in order to estimate only 
two states (X1Xo). 


The fourth-order continuous system is described by 


0 1 0 0 0 | [u(t)+w(t), . 
O O l 
x(t) = PAGE) Se | 
0 0 0 | 
-400 -320 -118 -19 soc 
A deterministic function u(t) = Of5 was used Gm the em 


ample. The random forcing function w(t) was assumed to 
have zero mean and variance equal to 0.01 square units, 
ieee ew eee CO 10 © 00 je 


The measurement equation is 
z(t) =;1 0 O Ojx(t) + v(t) , 


with the measurement noise v(t), considered to have zero 
mean and variance equal to 0.25 square units, i.e., 
v(t) = NiGOmOmeem2 5S )x 

Before applying the reduced-filter technique dis- 
cussed in Section III, A, the continuous process was dis- 
cretized using the ideas presented in Chapter II, so the 


discretized system could be described as 
x(k+1) = Ox(k) + Pfu(k) + wlk)] , 
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with measurement 
2(k) = {1 0 0 0] x(k) + v(k) , 


where 


ge" (2. 6a) 
was computed using the series expansion 


2 B 
dQ = I+ AT + si (AT) + ar(AT) + see . 
and 


Si 
me Gel es Cay Be (2.6b) 
O 


was computed by means of 


2 2 3 
eT Sap sr(AT) me. Gir - A = + AZ a =e 


A Monte Carlo simulation* was performed using 200 
replications of the trajectory, each trajectory being 


formed by 50 samples with sampling time T = 0.05 seconds. 


1 

An IBM System/360-67 computer was used for solving all 
the examples presented in this thesis. A copy of the 
main computer program for this example is presented in 
Appendix B. 


DS, 


This simulation was carried out in order to produce approx- 
imate ensemble values for the experimental mean and experi- 
mental variance of estimation error in xy and Xo, as a func- 
tion of time from t = 0, for both the complete, optimal fil- 
ter and the reduced filter. 

For the Monte Carlo simulation each state was assumed 
to have uncorrelated gaussian initial conditions with zero 


mean and the following variances: 


Ox, = 0.5 square units 
Z 
Tx, 
Z 
Ix, 
O*x 


2.0 square units 


2.0 square units 


fs 4.0 square units. 


In order to obtain an unbiased estimate from the 
Kalman filter, the initial values assigned to the filter 


state variables were 
&(0/-1) = E|x(0)| 


or 


2(0/-1) 
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and the initial value of P(0/-1l) was 


p(0/-1) = E [x(0) - %(0/-1) [x(0) - &(0/-1)] | 
[ 
if nae iT! 

= 4 | (0) . E[x(0) | | x(0 (0) = 1a Oy) 
\ 
= cov| x(0), F 
ong 
O55 
Cel 
Za) : 
eine Te 
LOC 
: 4.0 | 


The reduced filter was implemented to estimate only 
two states (x, and x,), and for this purpose partitioning 
was applied as follows. The ) and IT matrices (from 


Table III-1) 


OE he M0 0.0449 
-0.0066 Ol ESPARS) 


| 
| 
| 
(0.0 a= _.- . — ee... 
-0.3644 -0.2981 | 0.8871 0.0307 
| 
-12.2713 -10.1815 |-3.9181 0.3043 
T*(0.05) = [0.0001 0.0066 | 0.3644 12.2713], 
- } 


Sib 





were partitioned as shown to obtain 


0.9990 
d, = 


0.0066 


and 


Correspondingly, 


0.0449 
0.9946 


0.5 0 
Py ey ol) = | 
0) 20 


0 
¥)(0/-1) = ™ 
[9 J 
& 1 
ED ae 
and 
0.000 
Q) = 
| 0.000 


0.000 
0.000 


4 
| Gee taplertii-1) 


BZ 


The reduced filter was then implemented using the set 


of Eqs. (3.9a,b,c) repeated here for convenience, 
= eT) He! He Pe (k/k 1 yee & “(3.9a) 
Gj Cie) = Py i ae es | ated ac te (F298 


P,(k/k) = jI - G (k)Hy | Py (k/k-1) (3.9b) 


and 


ec / §,P, (k/k) 6, +Q. (3.9c) 


Attention Sou? Pe called ize the fact that the pre- 
sence of a ce agian txt ie cure Se in this problem 


changes slightly the form oe ve (3. r~ which can be re- 


is 
m3 i he 


written as 


X, (k/k) = 0, x, (k- -1/k- Ltr vu (ke 1)t G 14/200 
-H, ays yok L/k- La uw, (k- “»)) 


to reflect the effect of u(k) on the prediction step. 

The outputs of the computer simulation are shown in 
Tables II1I-2 through III-5, where numerical results versus 
Game ware tabulated for the gains, theoretical variange, 
experimental variance and experimental means of estimation 
error for both the complete and reduced filters. Figs. 
(3.3) through (3.6) show graphical representations com- 


paring complete and reduced filters performance. 


oe 


PHI MATRIX 


0.9999 
-0.0066 
-0. 3644 

22S 


GAMMA MATRIX 


0.0001 
0.0066 
0. 3644 
P22 13 


Q MATRIX 


0.0000 
0.0000 
0.0000 
0.0000 


0.0499 
0.9946 
-0.2981 
-10.1815 


0.0000 
0.0000 
0.0000 
0.0008 


Table III-l 


0.0012 
0.0480 
0.8871 


-3.9181 


0.0000 
0.0000 
0.0013 
0.0447 


0.0000 
0.0009 
0.0307 
0.3043 


0.0000 
0.0008 
0.0447 
1.5059 


VARIANCE OF FORCING NOISE = 0.01 
VARIANCE OF MEASUREMENT NOISE = 0.25 


INITIALIZATION - P(0/-1) 


0.5000 0.0 0.0 0.0 
0.0 2.0000 0.0 0.0 
0.0 0.0 2.0000 0.0 
0.0 0.0 0.0 4.0000 


INITIALIZATION - X(0/-1) 


0.0 
0.0 
0.0 
0.0 


T = 0.05 seconds. 
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Table [II-2 


GAINS FOR FULL AND REOUCED FILTERS 


G(2_,1) Gti3_,1) GI4,1) GRO1,1) GR2e1) 
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G (4x1) represents the matrix of gains for the optimal 


Pl biher. 


GR(2x1) represents the matrix of gains for the reduced 


filter. 
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Table III-3 


THEORETICAL VARTANCE OF ESTIMATION ERROR 
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gonal elements of the theoretical 


PK represents the dia 


filter. 


covariance matrix for the optimal 


duced filter. 


PR represents the diagonal elements of the theoretical 
covariance matrix for the re 
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Table [II-4 


EXPERIMENTAL VARIANCE OF ESTIMATION ERROR 
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VK represents the diagonal elements of the experimental 


covariance matrix for the optimal process. 


VR represents the diagonal elements of the experimental 


covariance matrix for the reduced process. 


37 

















Table III-5 


OF ESTIMATION ERROR 


EXPERIMENTAL MEANS 
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gh MEANE4 represent the mean of est 
error in the complete state vector by means of optimal 


MEANEI throu 
filtering. 


MEANEI1R and MEANE2R represent the mean of estimation 


error in states x, and x, by means of a reduced filter. 
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Variance of estimation error 
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Fig. 3.3 - Graphical plot of the variance of estima- 
tion error in x; vs. time. Curves a and b show the 
variance of estimation error, both theoretical and 
experimental, for the optimal estimation process. 
Curves c and d show the corresponding theoretical 
and experimental variances for the reduced case. 
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Fig. 3.4 - Graphical plot of the variance of est, 
error in x7 vs. time. Curves a and b show the vari- 
ance of estimation error, both theoretical and ex- 
perimental, for the optimal estimation process. 
Curves c and d show the corresponding theoretical 
and experimental variances for the reduced case. 
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Mean of estimation error 
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Fig. 3.5 - Graphical plot of the means of estimation 
error in xX} vs. time. Curve a shows the mean of 
estimation error for the optimal Kalman filter, and 
curve b shows that for the reduced filter. 
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Fig. 3.6 = Graphical plot of the means of estimation 

error in x> Vs. time. Gurvera shows the mean of 


estimation error for the optimal Kalman filter, and 
curve b shows that for the reduced filter. 
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IV. IMPLEMENTATION OF SUB-OPTIMAL LINEAR 
FRE leh oe hOReO och nar br NerilC TROCRSOES 
USING MATRIX PSEUDO-INVERS ION 


In the previous chapter a complete discussion of the 
implementation of a discrete, reduced model Kalman filter 
was presented, where the idea of partitioning the dynamics 
of the original plant was principal. In this chapter, a 
different approach to the same problem is presented, 
based upon the work of Meditch (Ref.3) for the continuous 
case. 

Applying the same idea of reduction of the original 
process, it is partitioned for this case into two or more 
reduced y; vector components (i = 1,2... M)y mith each 
new vector containing the same number of states. This 
new partitioning is accomplished by using matrix pseudo- 
inversion (see Appendix A). Once the original system has 
been broken down into Meereducedserder sub-optimal fil- 
ters, then the sub-optimal estimates for the original pro- 
cess are obtained by recombination of the outputs from all 
the sub-optimal filters. Section A presents the deriva- 
tion of the discrete sub-optimal filter formulation. 


These results are extended to Section B, where expressions 


1 
The order (n) of the original system is equal to M 


memes the order (r) of the reduced plants, that is 
n = Sie 
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for the sub-optimal estimates, theoretical sub-optimal 
estimation error and theoretical sub-optimal covariance of 
estimation error are derived. Finally in Section C an ex- 


ample is solved to illustrate the technigue. 


ae DERIVATION OF THE SUB-OPTIMAL FILTER EQUATIONS 

As in Chapter III-A it is assumed that the system 
dynamics are represented by Eq. (2.9), without consider- 
ing the presence of the deterministic part of the forcing 


vector, that is 

x(k+1) = Ox(k) + Dw(k) , 
with measurement 

gdom- Hx je vie 


where as before, w(k) and v(k) are random gaussian func- 
ELONS FILL yocdescrmibed byaboqs. (Z.ib) sebrougi ( 2amiege 

If, as discussed in Chapter I, it is desired to 
estimate only r states, where r<n, then the original vec- 
tor x(k) can be partitioned or transformed into two or 


more vectors by using a set of constant matrices 


ee 1. = | 2) eet 
a 
with A, having dimensions (r xn). Notice that each A. 
must have r rows and n columns. Then it is possible to 


define a new set of redticeamvect ors 
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CU as( Wane awl ie. . (4.1) 


with y; (k) pemaeicionitimeway (1m xa) vector. 

Now the concept of the Generalized Inverse? of a 
matrix (pseudo-inversion) is applied to the implementa- 
tion of the sub-optimal filters. A matrix A’ 2 is con- 
sidered to be the pseudo-inverse of a rectangular matrix 


A wf 
ANA =A, ie 


where in general A is a (s x t) matrix and he ie 4 Ge oe &) 
matrix. 

The selection of the set of constant matrices A, viele 
have as a necessary requirement for the solution of the 
filter implementation and the proper recovery of the sub- 


optimal estimates from the combination of the filter out- 


puts, that 
M 
) ATA, =I, (4.3) 
1=1 

I 


See Deutsch (Ref.10) and Penrose (Ref.8) for a complete 
discussion about matrix pseudo-inversion. 


2 


ats : 
A stands for the pseudo-inverse of A. 


This requirement will 


where I is the identity matrix. 
be clarified after the derivation which follows. 
If Eq. (4.1) is premultiplied by ae i=1,2 


and summation over all the i's is performed, then 


M M 

\ at ue ete 2h 

Ee = Ae 

i=1 i=] 
from which it is clearly seen, applying the requirement 
established by Eq. (4.3), that 

M 
+ 
Uns IS) (4.4) 
i=1 


This last equation provides the key for recovery of 


the original system state vector. 
Now if the original discrete system equation is pre- 


i=1,2...M, 


multiplied by A,, 
455) 


A,x(ktl) = A,Ox(k) + A, Dw(k) , 


and if Eqs. (4.1) and (4.4) are substituted am Eq. (4i5em 


the latter becomes 
M 
y; (k+1) ~ A. » Agy (k)+A, TD w(k), i=1,2 . M 
j=l 
(4.6) 


46 


Expression (4.6) represents a set of M reduced, dis- 
crete system equations, each of order r. 

The measurement process can be written assuming that 
all the measurements that were available for the complete 


filter are available too for the reduced filters. 
Zik) = Hxtk) + vee 
becomes after replacement eax k) by Eq. (4.4), 


M 
z(k) =H 2, Aty,(k) + w(k) (.47) 
re 


fe. (hal) willeheip to beeter visualeze the partitioned 
system. | 

In order to facilitate the following derivations, the 
optimal Kalman filter estimation Eq. (3.8) may be rewritten 


in a slightly different form, that is 


7 


R(k/k) = O2(k-1/k-1) + G(k)[2(k) - 2(k/k-1)] 
(4.8) 


where Z(k/k-1) represents the expected value of z(k) given 


all measurements up to time (k-1), or 


A(k/k-1) = Elz(k)| given z(k-1) 


4.7 


-uoTsSaaauT-opnasd xtTajzem Jo sueam Aq 'X smaysAs-qns pw o UT pauot3 


-Tj21ed maysks TeuT3TIO 9y} Jo uoTAeIUaSaidal MeAZeTP YOOTG - T°y “3TY 


fon 
aG 
—Z 
N| 
<{ 
| 
> 
Ne) 
re 
v 
am) 
| 
p> 
< 
(\ 


| M+ 
A) | Ss 4 OSE fy < 5 Aejyaq 3) ‘+ cS I 
t= 2 e% e% 
(7A) (T+) + /\ a (\ 


(4) 
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From Eq. (4.7), however, 


M 
2(k) =H) Ay yy (k) + v(k) (4.7) 
j=l 
Therefore, 
M 
2(k/k-1) = Hu) AG E y4(k)] given z(k-1) 
j=1 
or 
M 
B(k/k-1) =H) AT ,(k/k-1) (4.9) 
j=l 
Considering further that 
K(k/k-1) = QR(k-1/k-1) , Cero, 


where X(k/k-1) represents the expected value of x(k) given 
all measurements up to time (k-1l), and premultiplying Eq. 


ee 10) by A, 1,2... M, 
A,&(k/k-1) = A,O&(k-1/k-1) . (4.11) 
By extension of Eq. (4.1) 


A,&(k/k-1) = ¥,(k/k-1) , 
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and similarly by proper Use  ormugemC 
M 


van = i + aN 
x(k-1/k-1) = rae A. y,(k-1/k-1) 
= 


so Eq. (4.11) becomes 


M 
¥,(k/k-1) = A, ». AG Y,(k-1/k-1) . (4.1 
jz 


Basing the implementation of the reduced filters upon 
the dynamics of the set of systems described by Eq. (4.6), 
a set of M reduced-filter estimation equations may be 


stated by induction from Eq. (4.8), as 
9; (k/k) = 9,(k/k-1) + 6, (k)[2(k) - 2(k/k-1) | 
Replacing for y. (k/k-1) and Z2(k/k-1) Egs. (4.12) "and (4aa08 


gives 


M M 
9, (k/k)=a,9 ) AT (e-1/k-1) +6, ()| 2() -H \ F2,0/e-1)] 
j=l jel 
(4.13) 


W iG a elie eae 


This expression represents a set of M eA _ order coupled 


reduced filters, whose outputs ¥, (k/k) will be recombined 
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in a proper way in order to get the sub-optimal estimates 
x(k/k), where x (k/k) cman) |!) veckorm see Section 
IV-B). 

Referring to the optimal gain and theoretical vari- 
Mieewexpressions forwihe Kalman filter (Eqs. 3./7a,b,c), 
the equivalent M subsets of sub-optimal equations for each 
one of the reduced filters may be expressed, by definition, 


as 


7 T T 
a (ee P; (k/k-1)Hy (H,P; (k/k-1)H, +R (eas) 
where 
+ 
H, =H A, (see Eq. 4.7) 


eo) 


corresponds to the 1 subsystem contribution to the 


measurement vector; 

P,(k/k) =[I - 6, (k)H,|P.(k/k-1) , (4.14b) 
and 

P.(k+1/k) = 0,2, (k/k)O; +Q Ceili) 


where 


oo 
ee 
il 


+ 
A, oA, (see Eq. 4.6) 


A,DE[w(k)w (k)] TP "a; 


ro) 
il 


Tr 
A.QA; 


il 
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and 


It is important to notice that although the M estima- 
tion equations and the measurement equation are coupled 
for all the M reduced filters, the gain and covariance 


ee filter are truly uncoupled from Ge 


equations for the i 
other } 81,2 22. M filters, where sie: el 

Eqise (4.1 complete the implementation of the M sub- 
optimal filters. The filtering process will be visualized 
better by mneamomer sila cleo 

The initialization values for the sub-optimal filters 
are derived from the statistics used to initialize the com- 


plete optimal filter in the following way. Given the 


values for P(0/-1) and knowing that 
P,(0/-1) = Ely, (0)y; (0) 
iL (ti i 
and 


pinet 


P,(0/-1) = E A,x(0)x'(0)Al| 


A,E[x(0)x (0)] A; 


A.P(0/-1)A; (4.15) 


BV 


“(a/a)'S ‘laY[TJ paonpsazr yoses mozz sqnd3jno ayy 


ZutMous ‘SzaqTTJ Temtydo-qns pw 39uq Jo uot Ae{UasSvidar MeISeTP YS0TG - ZH *3TY 
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U (T-A/ADZ 


me ee ce “GE <n e ee  cmee eratin me a me, p> 


Es ee ee EE GSS ee ey 6 a ee) ee ee i ee te ca eR am 





By the same reasoning, 


y,(0/-1) = A,X(0/-1) (4.16) 


will be used as the initralization values for themset sone 


reduced estimation Faqs. (4.13). 


B. SUB=-OPTIMAL ESTIMATES, ESTIMATION ERROR AND SUB- 

OPTIMAL VARIANCE OF ESTIMATION ERROR 

In the preceding section, M reduced filters were im- 
plemented after breaking down the original process into M 
reduced (r x 1) vectors. In order to obtain the sub-optim- 
al estimates @t the original’ process, the outputs y, (k/k) 
from the reduced filters have to be recombined. This is 
easily accomplished by proper application of Eq. (4.4) to 


the filter's estimates, whereupon 


e 


M 
= = =“ fax 
x(k/k) = a Ay 94 (k/k) (4,17) 


will represent the sub-optimal estimates of the original 
process, with x(k/k) represa@m@imme an’ (mx IAvector. 

The remainder of this section is devoted to the de- 
velopment of a matrix difference equation for the itera- 
tive computation of the sub-optimal covariance of estima- 
Elon ero je 


If the error of sub-optimal estimation is defined as 


E(k) = x(k) - x(k/k) (4.18a) 
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then using Eq. (4.17) it becomes 


M 
E(k) = x(k) - ) AZ D,(k/k) (4.18b) 
j=l 
Or 
M 
e(kt+1l) = x(k+tl) - 2, AG 2, (kt1/k+1) eee ral sre ) 
cE 


eMPscicutnoneor Eqs..(4.83) and (4.7) into Eq. (4.18c) 


gives 


e(k+1l) = x(ktl) - 


M 
Nise 


M 
+r Cer aa. 
SUD EEC 


4 
Co 
Il 

— 


M M 
+6, (k+1)[H ie A y; (k+1)+v(k+1)-HO » at (k/te)] | 
j=l j=1 


M 

Replacing = Agy, (k+L) with x(k+l), and substituting 
ilaes 

x(k+1) = Ox(k) + Cw(k) from Eq. (3.1), then 


M M 
e(k+1) = Ox(k) + Tw(k) - 3S at] a8 » Aas (k/k) 
j=l i=! 
M 
+6, (k+1) [HOx(k) + HI w(k) + v(k+l) - HO 2, AFD (k/«)] | 
{= 


5)5) 


Expansion of this expression for the sub-optimal error 
leads to 


M M 
(k+l) = Ox(k) + Dw(k) - > aA.O we 


¥ 5 (k/k) 
jel 
M 
- ) ATG, (+1) HOx(k) - Yate, G.(k+1) HP w(k) 
el j=l 
oo us M 
ae 6, (k+1)v(k+1)+ Ag °C (+1) HO 2 
j= j 


Collecting terms, 


M 
e(kt1l) = 1G - 2 A, ve , (k+1)HO][ x(k) = Y- akp cure 


jal 


s eo 
AGG, (k+L)H | Pw(k)- ) AyGGc+1)vk+1). 


j=l 


Substitutienpiof Eqa(4518b) fopeeqhe cance 


M 
(k+l) = [I= AGG, (k+1)H!Oe(k)+ [1 my ee jDH\Dwk) 
‘lage j=1 
M 
- > AxG, (kt1)v(ktL). (4.19) 
j=l 
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Eq. (4.19) represents a theoretical form for the error 
of sub-optimal estimation which can be useful in establish- 
ing the steady-state estimation error for the sub-optimal 
Pilters. 

If the sub-optimal theoretical variance of estimation 


Seror ts defined as 
_ af= 7) 
S(k) = Ele(k)e (k)| (4.20) 


where 5(k) is an (nm x n)matrix, then substitution of Eq. 
(4.19) for e(k) gives 


f 


| M M 
S(k+1) = *{[p : 2, Aj; (+104) GEC» [I- 2, Ae oeeer ue 
iis I= 


aa r t 
a MATES j(erL (let) | | same vers | 


Carrying out the indicated matrix operations term by term 


yields 
‘ it 
- 7 + ~T ioe 
mBikt]) = E [1 2, Ag6y tH e(k)e (k)§ i ee 
im UF 


M M 
meets (k+1)H] TP w(k)w'(k)P '[1- mits A we (cH, 
l= j= 


e)iat G, a eee) G5 (+1) 45 "+ products 
j= —) | jal ) 


a 


Given the statistics of the problem, all the cross products 


will yield zero expectation. S(k+l) then becomes 


M M 
mn - T 
S(k+1) = [1 - ) ASG (k+l) H|OE fetkye '(k)) O° [ 1- » AG, k+])H | 
j=l j=l 
M M T 
+= ) AGG .(ctL)H]DE[w(k) w'QO] I [I- ) AGG ,(k+1)H | 
j=l j=l 


gE 


M 
‘ai + si 7. i ae ae 
* gaa wets (k+1)| aa (k+Da” , 
Je JT 


and finally 


M M 
. AE 
2 ; ani et oo ; 
S(k+ = i[ - Rot Gl eae Lye 5 (k i eG et lea) ee 
(k#L) = [T= ) AGG, (k+L)H OS@O'[T - ) AGG, (K+1)H I 
ae j=l 
M M T 
[ ee 1 eet 7 
Le Ge er ry eee let ) 
T= 2 A565( )H}Q{T - > AGG, (k+1)H | 
jel j=l 
M M 
- 3 A,G,(k+L)R | Ge(ktL)Ag (Lap 
j=1 j=l 


Eq. (4.21) represents an iterative expression for the 


computation of the theoretical sub-optimal covariance of 
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Po atiantcneseor lominrtialize the computation of S(kt1), 


the initial value of S(O) should be chosen as 


S(0) = E (0)e'(0) 


TH : 2 
Pea Gy al. 
Agy (O/-1) 


j=l 


ae CO) 


4 


E x(0) ~ K(0/-1) 7 


II 


PCy IL) 


The discrete matrix expression for the theoretical 
sub-optimal covariance of estimation error S(k), may be 
used as a reference to evaluate the performance of the re- 
ameecad set of filters, i.e., how close to the ‘trule values 


are the sub-optimal estimates. The solution for S(k) 


should be compared with the experimental variance of sub- 


optimal estimates. 


c. EXAMPLE 2 

The model used to investigate and to illustrate the 
technique presented in this chapter, was exactly the same 
fourth-order plant used in Example 1, subjected to the 
same deterministic and random forcing functions and with 
precisely the same measurements. The continuous system 


was discretized using the same sampling time T = 0.05 sec., 


52) 


and therefore the same @ and I’ matrices were obtained, 


arriving at the discretized system described by 

x(k+1) = Ox(k) + I } u(k) + wk), (2.79 
with measurement 

FAG 9 ital =0:10 OMA ®) 


The goal of this example, as in Example 1, was to 
estimate only two states of the fourth-order system 
(x7 and X) by means of a reduced-order filter. Apply- 
ing the technique discussed in Section A of this chapter, 
the original model (4 order) could be partitioned into 
two second-order models, but in order to make the compu- 
tational procedures much simpler, the original model was 
partitioned into four scalar models. This partitioning 
approach facilitated the selection of the A; matrices and 
the computation of their respective pseudo-inverses, A,. 

The solution to the partitioning problem was achiev- 
ed by trial and error. Four new vectors Yi» 1 sel 2, 3,4 


were defined by 
7; ae ES ep Dye 


The four A, (1 x 4) matrices selected were 
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It is important to remember here that all the A, matrices 
should have the same dimensions, and that the choice of the 
A's in this case renders the y,'s scalar quantities. 

The selection of the A. matrices was made by trial and 
error (in this case obvious and simple) in order to get the 
values of proper at's + that would fulfill the requirement 


established by Eq. (4.3), in this particular case 


The pseudo-inverse for each A; matrix was found using 


Bae (4.2) 


a a Ce2 


1 
Werther (Ref.9) presents a computation method that 
ylelds to matrix inversion or pseudo-inversion. 
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from whieh tones our Ay (4 x L) matrices were found to yo 


1 0 0 0) 
+ : + I + 0 + 

0. 0 le o 

0 | 0: 0 | a 


The partitioning of the original state vector is then 


summarized by 


gre Whe Omi VOC) = 24 (k) 
yy == 0 = Imre “OlxCky = x9 (k) 
Y3 = 0 O 1 0. x(k) = x3(k) 


y, = [0 0 O 1,x(k) = x,(k) 


with dynamics as specified by Eq. (4.6) 

This partitioning can be interpreted as an apparent 
canonization of the fourth-order model into four first- 
order models, represented each one by one state of the 
original system. 

Having selected the A. and Ay matrices, four 


coupled scalar filters were implemented in parallel using 


1 
t ¢an be easily verified, that the selection of the A. 
a A, matrices agrees with Eqs. (4.2) and (4.3). 
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Eqs. (4.13) and (4.14), but due to the presence of a deter- 
ministic forcing function, Eq. (4.13) was changed slightly 


as follows, 


4 
y.(k/k) = Ad) ay vy, (k- ae k=) EG. (ke) |2(k) 


ij 


4 


8 0 AG LY; (k/k-1) + Tyu(k-1), i=1,2,3,4 
jel 


Ca?) 
to show the effect of u(k) on the prediction step. 


The following dynamics were used for the computation 


of the filter gains and theoretical covariances. 


—= aE - ey 
d, = A, OA; ie IL chee 
‘fA a 
iL 0 0 Of ay Ae ae +p 1 
Vy = 22 0 
bee 0 
pee we 14g | Oy 
} 
fe lay , 


be 


where - 


and 


Also, 


where 


‘1 


Ds, 


and al are scalars. 


[fo D0 rere 
ee 

WD? 

[33 

ha 

Nat: i = eh 
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Similarly, 


and similarly, 


ip = $21 


J3™ 23) 
and 
Nee bay 


correspond to scalar dynamics of the reduced processes, 


-~ 


Note that the symbols (audit /9 we Ae Taea. 


whereas the symbols ane Mog. sere Mv yo recs CCE Cry: 
correspond to scalar elements of the dynamics of the full 
system. 


The Q matrix was partitioned as follows: 


— + Se 
Q, = A, QA; i= 1,2,8,4 
then 
| | 7 
i 9 0 off ayy qy Fi 
| | | 
7 ° 
q, = os on 
1 me le 
° | | 
94) 4419: 
. | 
a) a Saal j 
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and similarly 


19 3 0 
Ae - as ' 
and 
Sh Sulit, <i 
where Gis G2 +++ Gy, Go2 «ee etc., represent scalar quan- 
Eiplece 
The measurement process was specified by 
= . 
2(k) =H Agy 4 (k) + v(k) 
ce 
with the measurement matrix H = 1 O20 0}, and 


Vik). =: COM Oma O 2s) 
For gains and variance calculations for each filter, 
the measurement matrix H(l x 4) was partitioned as fol- 


lows: 





or 


h, = er 
hy = 0, 
h3 = QO, 
and 
hy = 0, 
with hy> hy ... being scalars. 


The initialization values for the reduced filters 
were computed using Eqs. (4.15) and (4.16), from which 
the initial guesses for the covariance of estimation 
error Pee) = 1) and expected value of the state vector 


y,(0/-1) were obtained as follows: 
Pa(0/-1) — A.P(O/-1)A, i = 1,2,3,4 


which gives 


i 
i 


1 0 0 0}! p,4(0/-1) fly 
Doo (0/-1) al iO: a 

p,(0/-1) = oe sa 
P33(0/-1) O 3 


t 


| b P4g6O/-L) | 0 


} 


6/7 


OL 


p, (0/-1) - p11 (0/-1) =O 550). 


ee See yal) = 2,0 | 


p3(0/-1) = p33(0/-1) = 2.0 , 
and 
Py, WO SIL) ae Dry 07-1) = a0 % 
where p, (0/-1) eet P17 60/-1) Sy ees. GORUerES Rates a: 


p,(0/-1), pp 07RD E oe: Cbce mer mceciibme nt hcaiind tla 
values of covariance for the reduced filters. 
For the initial values of the estimates it is easily 


seen that 


y, (0/-1) = &) (07-1) = 0.0; 


¥,(0/-1) = x,(0/-1) = 


| 
@, 
S 


y3(0/-1) = ,(0/-1) = 0.0 , 


and 


ll 
© 
© 


y,(0/-1) = x,(0/-1) 


will provide the filters with unbiased initial values. 
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a 


iene itour scalar, discrete, swh-optimal filters, in 
this example ecules only through the measurement process. 
were implemented and run in parallel in accordance with 
Eqs. (4.22), (4.14), and (4.7). At the end of each sam- 
pling interval the outputs from the four filters were com- 
bined by the expression 
4 
x(kK/k) =) AGP s(k/k) 


j=! 


e 


in order to recover the sub-optimal estimates of the 
original vector x(k), from which only the estimation of 
xX] and Xo was of particular interest. 

At the same time, the theoretical sub-optimal covari- 
PGCe Of EGtimation error Stige was computed using recur- 
Bieematrix Bq. (4.21), im order to evaluate the selection 
of the set of constant matrices Ay. This is possible to 
see when S(k) is compared with the experimental reduced 
covariance of estimation error. 

Again a Monte Carlo simulation* was performed with 


200 replications of the trajectory, each trajectory made 


iE 

This is not generally so. The uncoupling here is a con- 
sequence of the choice of the A; matrices. The yj state 
equations need not be uncoupled, but the associated gain 
and variance equations are, however, always uncoupled. 


Z 

As in Example 1, a computer System 360/67 was used for 
the simulation and solution of this problem. See Appendix 
C for corresponding computer program. 
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of 50 samples, with sampling intervals of 0.05 seconds. As 
a result of this simulation, ensemble values for the mean 
and experimental covariance of estimation error, as a func- 
tion of time, were obtained in Xy and xX» for the optimal 
and sub-optimal filters as a means of comparing their per- 
formance. 

Numerical results of this simulation are tabulated in 
Tables (IV-2) through (IV-6) and graphical representation 
of the results is shown in Figs, @G@mey) throucne(4.6)2 
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PHI MATRIX 


ORE, 
-0.0066 
-0.3644 

=12.2713 


GAMMA MATRIX 


0.0001 
0.0066 
0.3644 
eve? / 13 


Q MATRIX 


0.0000 
0.0000 
0.0000 
0.0000 


Oo 2 oo 


Table IV-1l 
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VARIANCE OF FORCING NOISE = 0.01 
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Table IV-2 


GAINS FOR FULL AND REDUCED FILTERS 
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represents the matrix of gains for the optimal 


filter. 


G (4x1) 


Ze 


Gl and G2 represent the gains used in the sub-optimal 
filters 1 and 


J? 
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VARTANCE OF ESTIMATION ERROR 
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PK represents the diagonal elements of the theoretical 


covariance matrix for the optimal filter. 
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perimental covariance matrix for the optimal 
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VK represents the diagonal elements of the ex- 
process. 
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SUBOPTIMAL VARTANCE OF ESTTUMATION ERROR 
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optimal theoretical covariance of estimation error, 


corresponding to the sub- 


SR represents two of the diagonal elements of the sub- 


optimal estimation of X} 


but it re- 
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y the same pattern as SR 


fers to the experimental covariance of estimation 


VR follows exactl 
error. 
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Table IV-6 


EXPERIMENTAL MEANS OF ESTIMATION ERROR 
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stimation 
f optimal 


MEANE] through MEANE4 represent the mean of e 


error in the complete state vector by means o 


filtering. 


MEANELR and MEANE2R represent the sub-optimal mean of 


1 and Xo by means of reduced 


estimation error in states x 


fiveers. 
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Variance of estimation error 





0.0 0.5 1.0 Ie) 20) ee) 
Time in seconds 


Fig. 4.3 - Graphical plot of the variance of estima- 
tion error in x} vs. time. Curves a and b show the 
experimental and theoretical variance of estimation 
error for the optimal process. Curves c and d show 
the sub-optimal experimental and theoretical co- 
variance of estimation error. 
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Variance of estimation error 
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Fig. 4.4 - Graphical plot of the variance of estima- 
tion error in x2 vs. time. Curves a and b show the 
experimental and theoretical variance of estimation 
error for the optimal process. Curves c and d show 
the sub-optimal case. 
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Mean of estimation error 
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Fig. 4.5 - Graphical plot of the mean of estimation 
error in x; vs. time. Curve a shows the mean of 

estimation error for the optimal estimation process. 
Curve b shows that for the sub-optimal case. 
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Meam of estimation error 
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Fig. 4.6 - Graphical plot of the means of estimation 
error in x9 vs. time. Curve a shows the mean of 
estimation error for the optimal estimation process. 
Curve b shows that for the sub-optimal case. 
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YE EMENT ALON, OR ASR EDUCED DISCRETE 
KALMAN FELTER BY USING THE DOMINANT 


ROOTS Or THEBOnRie INABS PROCESS 


The general idea of producing a reduced sub-optimal 
filter (Chapter III) or several sub-optimal filters of re- 
duced order (Chapter IV), is presented here using a dif- 


ferent approach. 


Ne THE DOMINANT-ROOTS APPROACH 
It 1s assumed as in the previous chapters, that the 


full system dynamics are represented by 

x(kt+1) = Ox(k) + P[u(k) + wlk)] , (2.9) 
with measurement 

z(k) = Hx(k) + v(k) , (2) 


and with the random processes involved, w(k) and v(k), 
fully described by Eqs. (2.11) through (2.16). 

It is clear that the order of the system is n and 
that the goal of the problem is to estimate only r states 
(r<n) of the process; therefore it will be necessary to 
produce a reduced model of the plant in order to use its 


dynamics in the implementation of the reduced filter. 
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It is common practice when handling classical auto- 
matic control problems to assume, or to require, that the 
high-order system be dominated by one pair of complex con- 
jugate roots. This will normally permit the use of a se- 
cond-order characteristic equation, and will simplify 
analysis and synthesis because the second-order system has 
been completely solved and all kind of curves, etc., are 


eee linear 


available. Applying the same idea to the n 
plant under consideration, but without the constraint of a 
pure second-order approximation, the dominant-roots repre- 
sentation will be determined mainly by the number r of 
states to be estimated and by the amount of performance de- 
gradation introduced when choosing the reduced model. 

It is noted that the number r of states to be estimat- 
ed will specify the minimum number of roots which can be 
selected as dominant, but an upper limit in the selection 
of the dominant roots would be determined by the designer's 
judgement, based upon the amount of degradation introduced 
by the reduced model. 

Due to the gualitative nature of the root-dominance 
criterion, it is very difficult to state a general proce- 
dure for the selection of the reduced-order model. There- 
fore a specific example will be used to illustrate the 


method. 


SZ 


Consider the fourth-order linear system with continu- 


ous dynamics used in Examples 1 and 2 and described by 


po 1 0 0 0 Vult)twlt)] 
. 0 0 1 0 0 
xc) = xe) ate 
0 0 0 : 
-400 -320 -118 -19) 400 | 
(oe) 
coeem 1. GC 
x(0) a QO ; 


and measurement 


z(t) = i o © 0] x(t) + v(t) 


The task of the problem is to estimate only two states 
of the original process (x, and Xo). 

When solving Example 1 (Chapter 3), the whole process 
(Eq. 5.1) was discretized and then the dynamics of the dis- 
crete system were partitioned according to the number of 
states to be estimated. Here, instead of discretizing the 
whole system, the order of the continuous process will be 
lowered by proper choice of the dominant roots, and then 


the reduced continuous system will be discretized. 
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The system described by Eq. (5.1) may be rewritten as 


x, (t) = Xo (t) ; (5.220) 
X(t) = X q(t) : ( 5, 2iy 
X(t) = xyalit) (5.2) 


and 
X(t) = -400x, (t)-320x5(t )-118x,(t)-19x, (t)+400 u(t )+tw(t). 
- (Seep 


For convenience a new variable J(t) = x, (t) is defined; 


EheneeGs wioeZa. bic) become 


Bc, 


= X5(t) 9 
7(t) = X(t) ’ 
o(t) oe K,(t) , 


and by direct substitution of these three equations into 


Eq. (5.2d) the Matter@becomes 


C(t) + 19 F(t) + 118 d(t) + 3204(t) + 4007(t) =400/ulc +w(t),, 


( Sep 


where 
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The dominant-root analysis of Eq. (5.3), which now fully 
represents the fourth-order system under consideration, 
will be made with the help of a pole-zero representation. 
This)wall require (asean easy approach) transformation of 
Eq. (5.3) into the s-domain, accomplished here by use of 
the Laplace transform. Then, considering all initial con- 


ditions equal to zero, Eq. (5.3) transforms into 
B(s)s*+19As) s°+118G(s) s“+3208(s)s+4008(s)=400 |U(s)+W(s) ]. 
which after proper factorization gives 


400lU(s)4+W(s) 


®(s) = (st10)(st5)(st2-j2)(st2+j2) ’ 


@)(s) 400 


U(s)+W(s) (st10)(st5)(st2-j2)(st2+j2) _ 


es): 


Eq. (5.4), may be interpreted as the transfer function of 
the fourth-order system relating ™(s) with the foxrerie 
functions (inputs). 

A pole-zero representation of Eq. (5.4) (Fig.5.1b), 
shows a pair of complex conjugate roots located at 
s = -2+j2, and two real roots located at s = -5 and s=-1l0 
respectively. If the dominance criterion is based upon the 


size of the real portion of the roots under consideration, 
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400 


(st+10)(st+5)(s*+4st8) na 











lad 
Ve juu 


Fig. 5.1 - a). Block diagram representation of the pro- 
cess described by Eq. (5.4). 

b). Pole-zero diagram for Eq. (5.4) indicating four roots 
at s = -5, s = -10, s = -24}2 respectively, and with 400 


considered to be the gain of the system. 


86 


then it is clear for this particular example, that the less 
dominant root will be at s = -10 (this will give rise to a 


gee and the most dominant roots will 


term of the form € 
be those at s = -2+}2, with the root at s = -5 producing an 
intermediate effect. 

As stated before, the purpose of this problem is to 
estimate only two states of the original system (xy and Ko). 
It therefore seems a good approach to select the roots at 
s = -2%}2 in order to produce a reduced second-order con- 
tinuous system, whose dynamics will be used in the imple- 
mentation of the reduced filter. It is important to note 
that the gain of the reduced system should be corrected in 
order to compensate for the disappearance of two of the 
roots of the original system. In this case the gain of the 


reduced system was made equal to the product of the magni- 


tudes of the dominant roots, that is 


gain reduction = |-2-32|/-2+j2| = 8. 


pee Fig. (5.2) 

Another important consideration should be based upon 
the fact that the performance degradation introduced, and 
the number m of inputs (forcing functions) to the original 
system, will strongly determine the order of the reduced 


model when using dominant roots. 
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® _(s) 








-jly 


Lae 


Fig. 5.2 - a). Block diagram representation of the 
reduced model. 
b). Pole-zero diagram showing only the dominant 


roots. Gain of the system has been adjusted to 8. 
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After the selection of the pair of complex conjugate 


roots at s = -2*j2, the reduced model may be expressed as 
©, (s) 2 8 _ 8 
U(s)+W(s) (st2-j2)(st+2+j2) s2+4st+8 


or transferring back to the differential equation forn, 
ae) ct 4Q(t) + 86(t) = 8fuCe)+w(t)] 


Recalling that E(t) = x, (t), the second-order reduced 


model becomes 


x4 (t) = Xo (t) 


Ko(t) = -8xy(t) -4xy(t) + 8[u(t)+w(t)), 


OF 


| 0 -1] °| [u(t )+w(e) | 
R(t) = | | Sg) . 
-4 


r ’ 
8 


which represents the reduced model from which the reduced 
filter will be implemented (see Fig. 5.2). 

No special derivation of the reduced-filter equations 
is necessary because the form of the equations will be ex- 
actly the same as those for the optimal filtering, and 


only the order will change. 


oh) 


The discrete reduced model is obtained from Eq. (5.5) 


using 
A_T 
Ss ; (5.6) 
and 
ALT -' -ALT 
r=€ ap € Y aye, (5.7) 
r ie 
O 


with Gs ye =) 'shnla BL(r x m) representing the reduced sys- 
tem and reduced distribution matrices respectively for the 
continuous case, and d(x x Yr) and Ie a x ™) being the 
state transition matrix and distribution matrix for the 
reduced discrete case. Then the discretized reduced pro- 


cess can be described as 


ee Koel) ec Tf u(k )+w(k)] 


with measurement 


z(k) = Hx(k) + v(k) 


As in the partitioned method presented in Chapter 3, 
all the measured outputs are considered to be included in 
the reduced model, and furthermore the H matrix is assumed 


to be partitioned as 


with H, = 0, and H, (p omic ye 


20 


Considering the problem of initialization for the dis- 
crete reduced covariance equations, it is solved exactly as 
in Chapter 3, by partitioning the P(0/-1) matrix from the 


filter, as 


[Py (0/-1) | Py2(0/-1) | 


PCO/slL) =|) eS Se 


and choosing 


Ai) Po Oy) 


with Po = (r x r) and ea Oy all) Geax). the latter 
representing the initial value of covariance for the re- 
duced filter. 


Hor the teduced filter initial states, 
ZN _ A 7 
R(0/-1) = @,(0/-1) 


where 


x, (0/-1) 
Ro/-1) = |-- --- 
X5(0/-1) 
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sa EXAMPLE 3 
The same example used in the previous methods is em- 
ployed here to illustrate the method of dominant roots. 
The same initialization values, same statistics for 
the random sequences, and same measurements as those used 
in Example 1, were used im this exapple for both the fui 
and reduced systems. The corresponding values for the re- 


duced filter are listed here Eor Convenwence 


"0 | 10.5 0 | 
| 
se (Q0)/ ly) = ‘ P_(O7PL) = j 
op 0120 
H= (1 O| , 


and 


vo mcmniGO. 0, 0.25) . 


wik) 2 NCORG 0201) 


A significant difference between the two examples is the 
determination of the reduced transition and distribution 
matrices which in this case were computed as indicated by 


Bas... (.5.6 and (S77) wel 


a 


Due to the similarity between the partitioning method and 
the dominant-roots approach for this particular problen, 
the same computer program was used for the Monte Carlo simula- 
tion and solution of Examples 1 and 3, with the only change 
introduced regarding the computation of the dynamics for 
the reduced model. 

Table (V-1) shows the initialization values for the 
reduced model, as well as the values for the O and I’. ma- 
trices. Tables (V-2) through (V-5) present numerical re- 
sults of the simulation, and a graphical representation of 
the simulation is shown in Figs. (5.3) through (5.6). 

Fig. (5.7) shows a complete block diagram of the re- 


duced filtering process for this particular example. 
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PHI MATRIX 


On 5999 
-0.0066 
-0.3644 

-12 22713 


GAMMA MATRIX 


0.0001 
0.0066 
0.3644 
12.27% 


REDUCED PHI 


0.9906 
-0.3613 


Table V-1l 


0.0499 
0.9946 
-0.2981 
-10.1815 


MATRIX 


0.0452 
0.8100 


REDUCED GAMMA MATRIX 


O70093 
0.3613 


Q MATRIX (REDUCED) 


0.0000 
0.0000 


VARIANCE OF FORCING NOISE 
VARIANCE OF MEASUREMENT NOISE = 0.25 


0.0000 
0.0000 


REDUCED P(0/-1) 


0.5000 
0.0 


INITIALIZATION - X(0/-1) 


0.0 
0.0 
0.0 
O.0 


0.0 
2.0000 


0.0012 
0.0480 
0.8871 
-3.9181 


= 0.01 
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0.0000 
0.0009 
020307 
0.3043 


Table V-2 


GAINS FOR FULL AND REDUCED FILTERS 


G(2,1) G(3,1) G(4,1) GR(1,1) GR251) 


G(1l,1) 
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Coeoeoeeoeroer eee eee @ 
C2222 0° 
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G (4xl) represents the matrix of gains for the optimal 


filter. 


GR(2x1) represents the matrix of gains for the reduced 


filter. 
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Table V-3 


THEORETICAL VARIANCE OF ESTIMATION ERROR 


PK(2,2) PK(3,3) PK(4_4) PR( 1,1) PR(292) 


PK(1,1) 
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PK represents the diagonal elements of the theoretical 


covariance matrix for the optimal filter. 


PR represents the diagonal elements of the theoretical 


covariance matrix for the reduced filter. 


96 


Table V-4 


EXPERIMENTAL VARIANCE OF ESTIMATION ERROR 


VK(292) VK (3,3) VK(4,4) VRO1s1) VRI292) 


VK(1,1) 
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VK represents the diagonal elements of the experimental 


covariance matrix for the optimal process. 


VR represents the diagonal elements of the experimental 


covariance matrix for the reduced process. 
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Table V-5 


EXPERIMENTAL MEANS OF ESTIMATION ERROR 
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error in the complete state vector by means of optimal 


MEANE1] through MEANE4 represent the mean of estimation 
filtering. 


MEANE1R and MEANE2R represent the mean of estimation 


and x» by means of a reduced filter. 
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Fig. 5.3 - Graphical plot of the variance of estima- 
tion error in x; vs. time. Curves a and b show the 
theoretical variance of estimation error for the 
reduced and full filters respectively. Curves c and 
d show the experimental variance for the full and 
reduced processes respectively. 
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Fig. 5.4 - Graphical plot of the variance of estima- 
tion error in x37 vs. time. Curves a and b show the 
theoretical variance of estimation error for the 
reduced and full filters respectively. Curves c and 
d show the experimental variance for the full and 
reduced processes respectively. 
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Fig. 3.5 - Graphical plot of the means of estimation 
error in xX} vs. time. Curve a shows the mean of 
estimation error for the optimal Kalman filter, and 
curve b shows that for the reduced filter. 
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Petes e). bums Graphical plot of the means of estimation 
error in x9 vs. time. Curve a shows the mean of 


estimation error for the optimal Kalman filter, and 
curve b shows that for the reduced filter. 
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VI. DISCUSS TONSORS RES Uris 


The three different approaches presented in this 
thesis have as a common objective a decrease in the compu- 
tational requirements involved in handling systems of high 
dimension, based upon the fact that only a preselected re- 
duced number (r) of states is of significant importance in 
the estimation process. 

For the particular example chosen as an illustration 
for the three methods, the partitioning idea discussed in 
Chapter III, although it is the easiest to apply, appears 
to give the least attractive results in the reduced esti- 
mation of states xy and Xo where a small bias is produced. 
This seems to be a direct consequence of the degradation 
introduced when 0.4 (which in this example has relative 
large elements as compared with 0,) was not considered to 
be a part of the dynamics of the primary system. 

In the pseudo-inversion approach, implementation of 
four scalar filters, plus the selection of a proper set of 
A, matrices and its pseudo-inverses Ay was necessary, but 
Since it takes less computing capacity to manipulate four 
scalar filters than one (4 x 4) #ilter, the estimation 


process was clearly simplified. 
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The selection of the partitioning matrices A, is arbi- 
brat hon eanyepant LcUleipmoecn, ebinerefore for this ex- 


ample another selection was made as 


Ay = | 1 -l -l , 
‘ ] 
a oleelaeesl 
ie 
A,= 1 -1 -1 1], 
and 
A,={l -2 1 -1| 
with 
owas /0.25' PO. 25, 0.25] 
| : 
Oe 5 On2s9 On? 5 | ZnS 
he = ee | ; ie re | 9 At =a | ; AY i. 
ie -0.25 2 bee 3 -0.25 4 0.25, 
| : : a 
0 asy 0.25) | 0.25] |-0.25) 


This particular choice, however, leads to four fil- 
ters for which the sub-optimal variance of estimation er- 
ror, both theoretical and experimental, increases linearly 
with time, which is undesireable. Therefore this parti- 
cular partitioning was not used. 

The plant for all three examples was selected with a 
pair of clearly dominant complex roots in order to obtain 
a good second-order approximation of the complete process 
from which the reduced filter could be implemented using 


the method presented in Chapter V. 


HO 


APPENDIX A 


MATRIX PSEUDO INVERSION 
A review of matrix pseudo-inversion is presented in: 
this appendix as an extension to the idea of pseudo-inver- 


sion introduced in Chapter IV. 


, a : : 
A matrix A is said to be the pseudo-inverse of a rec- 


tangular matrix A if 


e e r o @ — L 
where if A is a square non-singular matrix, i.e., A ex- 


ists, then 


1s 4 trivial s@lkitione since A Neale =< ie 


Another important definition may be introduced using 


a set of four matrix equations 


AAT A=A , (A-2a) 
Tien aie (A=2ip 
T 
A at|" = aat, (A-2c) 
and 
T 
[Ata] = Ata (A-2d) 


106 


en — 


This definition implies that the inverse of a matrix A is 
established as the unique pseudo-inverse if and only if it 
satisfies the set of Eqs. (A-2). Furthermore, as a conse- 
quence of the latter definition, the pseudo-inverse among 


other properties has the following: 


++ 


A = A 

alt = nen 

4s -1 

A = A ~, when A is a non-singular matrix, 
and 

7 5 Oe 


The pseudo-inverse AY of a matrix A may be shown to 
be unique as follows. 
A : 
Let A. be a set of pseudo-inverse of A. where 


: + 
i1=1,2 ...M, and assume that the A. are chosen so that 


M 
ale - 
» ata, = I. (A-3) 


+ 


Further assume that another set of pseudo-inverse ve * A; 


exists for each Ass and that this new set of e also sat- 


isoties 


M 
> ta, =. (A-4) 
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Then it follows by proper use of Eqs. (A-3) and (A-4) that 


and since the choice of A, is arbitrary, it follows that 
Y. = A. with i = 1,2). 


and the uniqueness property is established. 
For more complete treatment of this subject see 


Penrose (Ref.8), Werther (Ref.9) and Deutsch (Ref.10). 
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